## Load Packages ##
using Pkg
using Plots
using QuadGK
using LaTeXStrings

## Set paths ##
#path_out = "C:/Users/f006n8h/Dropbox/00_Transparency_project/04_output"
path_out = "/Users/F006N8H/Dropbox/00_Transparency_project/04_output"


## Set plot style ##
theme(:default) # Reset to default theme to start with a white background
plot!(background_color=:white) # Ensure the plot background is white


## Set functions ##

function z(mu, λ, z_max)
    if mu < λ
        0.0
    else
        z_max * (mu - λ) / (1.0 - λ)
    end
end

function integral_part(p_bar, v, mu, N)
    function integrand(p)
        ((v - p) / p)^(1 / (N - 1))
    end
    integral, _ = quadgk(integrand, p_bar, v)
    integral
end

function E_p(mu, N, v)
    p_bar = v / ((mu * N / (1 - mu)) + 1)
    ((1 - mu) / (N * mu))^(1 / (N - 1)) * integral_part(p_bar, v, mu, N) + p_bar
end

function E_p_min(mu, E_p_value)
    (1 - mu) / mu * (v - E_p_value)
end

function numerical_derivative(f, mu, h=1e-5)
    (f(mu + h) - f(mu - h)) / (2*h)
end


## Figure 1, (a) and (b) ##

# Parameters for Fig1a
N = 2
v = 1.0
λ = 0.2
z_max = 0.5
epsilon = 0.000001
mu_values = epsilon:0.000001:1-epsilon

# Calculating values for Fig1a plots
z_values = [z(mu, λ, z_max) for mu in mu_values]
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values = E_p_values - E_p_min_values

# Find the μ value that maximizes VOI
max_VOI_index = argmax(VOI_values)
max_VOI_mu = mu_values[max_VOI_index]
max_VOI_value = maximum(VOI_values) * 1.05 # Adjusted for better annotation placement

# Plotting Fig1a
Fig1a = plot(mu_values, z_values, label="Access Cost z", color=:blue, linewidth=2, xlabel="μ", ylabel="", legend=:topleft, size=(800, 600), framestyle=:origin, legendfontsize=14, tickfontsize=14)
plot!(mu_values, VOI_values, label="Value of Information (VOI)", color=:red, linewidth=2, linestyle=:dash)

# Vertical line for max_VOI_mu
max_y_for_line = maximum(VOI_values) # Adjust as necessary
plot!([max_VOI_mu, max_VOI_mu], [0.0, max_y_for_line], color=:black, linestyle=:solid, linewidth=2, label="")

# Adjusting annotation to appear on the x-axis
annotate!(Fig1a, max_VOI_mu, max_y_for_line , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu, digits=2))"), 14, :center, :bottom))

# Show plot
display(Fig1a)

# Export to PDF
savefig(Fig1a, join([path_out, "Fig1a.pdf"], "/"))




# Parameters for Fig1b
N = 5
# v, λ, and z_max remain the same as in Fig1a

# Calculating values for Fig1b plots
z_values = [z(mu, λ, z_max) for mu in mu_values]
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values = E_p_values - E_p_min_values

# Find the μ value that maximizes VOI
max_VOI_index = argmax(VOI_values)
max_VOI_mu = mu_values[max_VOI_index]
max_VOI_value = maximum(VOI_values) * 1.05 # Adjusted for better annotation placement

# Plotting Fig1b
Fig1b = plot(mu_values, z_values, label="Access Cost z", color=:blue, linewidth=2, xlabel="μ", ylabel="", legend=:topleft, size=(800, 600), framestyle=:origin, legendfontsize=14, tickfontsize=14)
plot!(mu_values, VOI_values, label="Value of Information (VOI)", color=:red, linewidth=2, linestyle=:dash)

# Vertical line for max_VOI_mu
max_y_for_line = maximum(VOI_values) # Adjust as necessary
plot!([max_VOI_mu, max_VOI_mu], [0.0, max_y_for_line], color=:black, linestyle=:solid, linewidth=2, label="")

# Adjusting annotation to appear on the x-axis
annotate!(Fig1b, max_VOI_mu, max_y_for_line , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu, digits=2))"), 14, :center, :bottom))

# Show plot
display(Fig1b)

# Export to PDF
savefig(Fig1b, join([path_out, "Fig1b.pdf"], "/"))





## Appendix Figure 1, (a), (b), (c), (d) ##

# Parameters for App Fig1a
N = 2
λ = 0.2
z_max = 0.5
epsilon = 0.000001
mu_values = epsilon:0.000001:1-epsilon

# Calculating values for App.Fig1a plots
z_values = [z(mu, λ, z_max) for mu in mu_values]

v = 1.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_1 = E_p_values - E_p_min_values

v = 2.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_2 = E_p_values - E_p_min_values

v = 5.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_5 = E_p_values - E_p_min_values


# Find the μ value that maximizes VOI
max_VOI_index_1 = argmax(VOI_values_1)
max_VOI_mu_1 = mu_values[max_VOI_index_1]
max_VOI_value_1 = maximum(VOI_values_1) * 1.05 # Adjusted for better annotation placement
max_VOI_index_2 = argmax(VOI_values_2)
max_VOI_mu_2 = mu_values[max_VOI_index_2]
max_VOI_value_2 = maximum(VOI_values_2) * 1.05 # Adjusted for better annotation placement
max_VOI_index_5 = argmax(VOI_values_5)
max_VOI_mu_5 = mu_values[max_VOI_index_5]
max_VOI_value_5 = maximum(VOI_values_5) * 1.05 # Adjusted for better annotation placement

# Plotting Appendix Fig1a
AppFig1a = plot(mu_values, z_values, label="Access Cost z", color=:blue, linewidth=2, xlabel="μ", ylabel="", legend=:topleft, size=(800, 600), framestyle=:origin, legendfontsize=14, tickfontsize=14, ylims=(0,3))
plot!(mu_values, VOI_values_1, label="VOI v = 1", color=:red, linewidth=2, linestyle=:dash)
plot!(mu_values, VOI_values_2, label="VOI v = 2", color=:green, linewidth=2, linestyle=:dot)
plot!(mu_values, VOI_values_5, label="VOI v = 5", color=:purple, linewidth=2, linestyle=:dashdot)

# Vertical line for max_VOI_mu
max_y_for_line_1 = maximum(VOI_values_1) # Adjust as necessary
max_y_for_line_2 = maximum(VOI_values_2) # Adjust as necessary
max_y_for_line_5 = maximum(VOI_values_5) # Adjust as necessary
#plot!([max_VOI_mu, max_VOI_mu], [0.0, max_y_for_line], color=:black, linestyle=:solid, linewidth=2, label="")

# Adjusting annotation to appear on the x-axis
annotate!(AppFig1a, max_VOI_mu_1, max_y_for_line_1 , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu_1, digits=2))"), 14, :center, :bottom))
annotate!(AppFig1a, max_VOI_mu_2, max_y_for_line_2 , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu_2, digits=2))"), 14, :center, :bottom))
annotate!(AppFig1a, max_VOI_mu_5, max_y_for_line_5 , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu_5, digits=2))"), 14, :center, :bottom))

# Show plot
display(AppFig1a)

# Export to PDF
savefig(AppFig1a, join([path_out, "AppFig1a.pdf"], "/"))



# Parameters for App Fig1b
N = 3
λ = 0.2
z_max = 0.5
epsilon = 0.000001
mu_values = epsilon:0.000001:1-epsilon

# Calculating values for App.Fig1a plots
z_values = [z(mu, λ, z_max) for mu in mu_values]

v = 1.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_1 = E_p_values - E_p_min_values

v = 2.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_2 = E_p_values - E_p_min_values

v = 5.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_5 = E_p_values - E_p_min_values


# Find the μ value that maximizes VOI
max_VOI_index_1 = argmax(VOI_values_1)
max_VOI_mu_1 = mu_values[max_VOI_index_1]
max_VOI_value_1 = maximum(VOI_values_1) * 1.05 # Adjusted for better annotation placement
max_VOI_index_2 = argmax(VOI_values_2)
max_VOI_mu_2 = mu_values[max_VOI_index_2]
max_VOI_value_2 = maximum(VOI_values_2) * 1.05 # Adjusted for better annotation placement
max_VOI_index_5 = argmax(VOI_values_5)
max_VOI_mu_5 = mu_values[max_VOI_index_5]
max_VOI_value_5 = maximum(VOI_values_5) * 1.05 # Adjusted for better annotation placement

# Plotting Appendix Fig1b
AppFig1b = plot(mu_values, z_values, label="Access Cost z", color=:blue, linewidth=2, xlabel="μ", ylabel="", legend=:topleft, size=(800, 600), framestyle=:origin, legendfontsize=14, tickfontsize=14, ylims=(0,3))
plot!(mu_values, VOI_values_1, label="VOI v = 1", color=:red, linewidth=2, linestyle=:dash)
plot!(mu_values, VOI_values_2, label="VOI v = 2", color=:green, linewidth=2, linestyle=:dot)
plot!(mu_values, VOI_values_5, label="VOI v = 5", color=:purple, linewidth=2, linestyle=:dashdot)

# Vertical line for max_VOI_mu
max_y_for_line_1 = maximum(VOI_values_1) # Adjust as necessary
max_y_for_line_2 = maximum(VOI_values_2) # Adjust as necessary
max_y_for_line_5 = maximum(VOI_values_5) # Adjust as necessary
#plot!([max_VOI_mu, max_VOI_mu], [0.0, max_y_for_line], color=:black, linestyle=:solid, linewidth=2, label="")

# Adjusting annotation to appear on the x-axis
annotate!(AppFig1b, max_VOI_mu_1, max_y_for_line_1 , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu_1, digits=2))"), 14, :center, :bottom))
annotate!(AppFig1b, max_VOI_mu_2, max_y_for_line_2 , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu_2, digits=2))"), 14, :center, :bottom))
annotate!(AppFig1b, max_VOI_mu_5, max_y_for_line_5 , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu_5, digits=2))"), 14, :center, :bottom))

# Show plot
display(AppFig1b)

# Export to PDF
savefig(AppFig1b, join([path_out, "AppFig1b.pdf"], "/"))



# Parameters for App Fig1c
N = 5
λ = 0.2
z_max = 0.5
epsilon = 0.000001
mu_values = epsilon:0.000001:1-epsilon

# Calculating values for App.Fig1a plots
z_values = [z(mu, λ, z_max) for mu in mu_values]

v = 1.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_1 = E_p_values - E_p_min_values

v = 2.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_2 = E_p_values - E_p_min_values

v = 5.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_5 = E_p_values - E_p_min_values


# Find the μ value that maximizes VOI
max_VOI_index_1 = argmax(VOI_values_1)
max_VOI_mu_1 = mu_values[max_VOI_index_1]
max_VOI_value_1 = maximum(VOI_values_1) * 1.05 # Adjusted for better annotation placement
max_VOI_index_2 = argmax(VOI_values_2)
max_VOI_mu_2 = mu_values[max_VOI_index_2]
max_VOI_value_2 = maximum(VOI_values_2) * 1.05 # Adjusted for better annotation placement
max_VOI_index_5 = argmax(VOI_values_5)
max_VOI_mu_5 = mu_values[max_VOI_index_5]
max_VOI_value_5 = maximum(VOI_values_5) * 1.05 # Adjusted for better annotation placement

# Plotting Appendix Fig1c
AppFig1c = plot(mu_values, z_values, label="Access Cost z", color=:blue, linewidth=2, xlabel="μ", ylabel="", legend=:topleft, size=(800, 600), framestyle=:origin, legendfontsize=14, tickfontsize=14, ylims=(0,3))
plot!(mu_values, VOI_values_1, label="VOI v = 1", color=:red, linewidth=2, linestyle=:dash)
plot!(mu_values, VOI_values_2, label="VOI v = 2", color=:green, linewidth=2, linestyle=:dot)
plot!(mu_values, VOI_values_5, label="VOI v = 5", color=:purple, linewidth=2, linestyle=:dashdot)

# Vertical line for max_VOI_mu
max_y_for_line_1 = maximum(VOI_values_1) # Adjust as necessary
max_y_for_line_2 = maximum(VOI_values_2) # Adjust as necessary
max_y_for_line_5 = maximum(VOI_values_5) # Adjust as necessary
#plot!([max_VOI_mu, max_VOI_mu], [0.0, max_y_for_line], color=:black, linestyle=:solid, linewidth=2, label="")

# Adjusting annotation to appear on the x-axis
annotate!(AppFig1c, max_VOI_mu_1, max_y_for_line_1 , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu_1, digits=2))"), 14, :center, :bottom))
annotate!(AppFig1c, max_VOI_mu_2, max_y_for_line_2 , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu_2, digits=2))"), 14, :center, :bottom))
annotate!(AppFig1c, max_VOI_mu_5, max_y_for_line_5 , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu_5, digits=2))"), 14, :center, :bottom))

# Show plot
display(AppFig1c)

# Export to PDF
savefig(AppFig1c, join([path_out, "AppFig1c.pdf"], "/"))



# Parameters for App Fig1d
N = 10
λ = 0.2
z_max = 0.5
epsilon = 0.000001
mu_values = epsilon:0.000001:1-epsilon

# Calculating values for App.Fig1d plots
z_values = [z(mu, λ, z_max) for mu in mu_values]

v = 1.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_1 = E_p_values - E_p_min_values

v = 2.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_2 = E_p_values - E_p_min_values

v = 5.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_5 = E_p_values - E_p_min_values


# Find the μ value that maximizes VOI
max_VOI_index_1 = argmax(VOI_values_1)
max_VOI_mu_1 = mu_values[max_VOI_index_1]
max_VOI_value_1 = maximum(VOI_values_1) * 1.05 # Adjusted for better annotation placement
max_VOI_index_2 = argmax(VOI_values_2)
max_VOI_mu_2 = mu_values[max_VOI_index_2]
max_VOI_value_2 = maximum(VOI_values_2) * 1.05 # Adjusted for better annotation placement
max_VOI_index_5 = argmax(VOI_values_5)
max_VOI_mu_5 = mu_values[max_VOI_index_5]
max_VOI_value_5 = maximum(VOI_values_5) * 1.05 # Adjusted for better annotation placement

# Plotting Appendix Fig1d
AppFig1d = plot(mu_values, z_values, label="Access Cost z", color=:blue, linewidth=2, xlabel="μ", ylabel="", legend=:topleft, size=(800, 600), framestyle=:origin, legendfontsize=14, tickfontsize=14, ylims=(0,3))
plot!(mu_values, VOI_values_1, label="VOI v = 1", color=:red, linewidth=2, linestyle=:dash)
plot!(mu_values, VOI_values_2, label="VOI v = 2", color=:green, linewidth=2, linestyle=:dot)
plot!(mu_values, VOI_values_5, label="VOI v = 5", color=:purple, linewidth=2, linestyle=:dashdot)

# Vertical line for max_VOI_mu
max_y_for_line_1 = maximum(VOI_values_1) # Adjust as necessary
max_y_for_line_2 = maximum(VOI_values_2) # Adjust as necessary
max_y_for_line_5 = maximum(VOI_values_5) # Adjust as necessary
#plot!([max_VOI_mu, max_VOI_mu], [0.0, max_y_for_line], color=:black, linestyle=:solid, linewidth=2, label="")

# Adjusting annotation to appear on the x-axis
annotate!(AppFig1d, max_VOI_mu_1, max_y_for_line_1 , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu_1, digits=2))"), 14, :center, :bottom))
annotate!(AppFig1d, max_VOI_mu_2, max_y_for_line_2 , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu_2, digits=2))"), 14, :center, :bottom))
annotate!(AppFig1d, max_VOI_mu_5, max_y_for_line_5 , text(latexstring("\\hat{\\mu} = $(round(max_VOI_mu_5, digits=2))"), 14, :center, :bottom))

# Show plot
display(AppFig1d)

# Export to PDF
savefig(AppFig1d, join([path_out, "AppFig1d.pdf"], "/"))







## Figure 2, (a) and (b) ##

# Parameters for Fig2a
N = 2
v = 1.0
λ = 0.2
z_max = 0.5
epsilon = 0.05
mu_values = epsilon:0.01:1-epsilon

# Calculating values for Fig1a plots
z_values = [z(mu, λ, z_max) for mu in mu_values]
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values = E_p_values - E_p_min_values
E_p_derivative_values = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]

# Find the μ value that maximizes VOI
max_VOI_index = argmax(VOI_values)
max_VOI_mu = mu_values[max_VOI_index]
max_VOI_value = maximum(VOI_values) * 1.05 # Adjusted for better annotation placement

# Find the μ value that maximizes E_p_derivative
max_Epd_index = argmax(E_p_derivative_values)
max_Epd_mu = mu_values[max_Epd_index]
max_Epd_value = maximum(E_p_derivative_values) * 1.05 # Adjusted for better annotation placement

# Find the μ value that maximizes E_p_min_derivative
max_Epmd_index = argmax(E_p_min_derivative_values)
max_Epmd_mu = mu_values[max_Epmd_index]
max_Epmd_value = maximum(E_p_min_derivative_values) * 1.05 # Adjusted for better annotation placement

# Plotting Fig2a
Fig2a = plot(mu_values, VOI_values, label="Value of Information (VOI)", color=:red, linewidth=2, linestyle=:solid, xlabel="", ylabel="", legend=:bottom, size=(800, 600), framestyle=:origin, legendfontsize=14, tickfontsize=14)
plot!(mu_values, E_p_derivative_values, label="dE[p]/dμ", color=:blue, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_min_derivative_values, label="dE[p_min]/dμ", color=:green, linewidth=2, linestyle=:dashdot)

# Vertical line for max_VOI_mu
max_y_for_line_voi = maximum(VOI_values) # Adjust as necessary
max_y_for_line_Epd = maximum(E_p_derivative_values) # Adjust as necessary
max_y_for_line_Epmd = maximum(E_p_min_derivative_values) # Adjust as necessary
#plot!([max_VOI_mu, max_VOI_mu], [-1.0, max_y_for_line], color=:black, linestyle=:solid, linewidth=2, label="")

# Adjusting annotation to appear on the x-axis
annotate!(Fig2a, max_VOI_mu, max_y_for_line_voi , text(latexstring("$(round(max_VOI_mu, digits=2))"), 14, :center, :bottom))
annotate!(Fig2a, max_Epd_mu, max_y_for_line_Epd , text(latexstring("$(round(max_Epd_mu, digits=2))"), 14, :center, :bottom))
annotate!(Fig2a, max_Epmd_mu, max_y_for_line_Epmd , text(latexstring("$(round(max_Epmd_mu, digits=2))"), 14, :center, :bottom))
annotate!(Fig2a, 0.95, -0.03, text("μ", :right, 14))


# Show plot
display(Fig2a)

# Export to PDF
savefig(Fig2a, join([path_out, "Fig2a.pdf"], "/"))


# Parameters for Fig2b
N = 5

# Calculating values for Fig2b plots
z_values = [z(mu, λ, z_max) for mu in mu_values]
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values = E_p_values - E_p_min_values
E_p_derivative_values = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]

# Find the μ value that maximizes VOI
max_VOI_index = argmax(VOI_values)
max_VOI_mu = mu_values[max_VOI_index]
max_VOI_value = maximum(VOI_values) * 1.05 # Adjusted for better annotation placement

# Find the μ value that maximizes E_p_derivative
max_Epd_index = argmax(E_p_derivative_values)
max_Epd_mu = mu_values[max_Epd_index]
max_Epd_value = maximum(E_p_derivative_values) * 1.05 # Adjusted for better annotation placement

# Find the μ value that maximizes E_p_min_derivative
max_Epmd_index = argmax(E_p_min_derivative_values)
max_Epmd_mu = mu_values[max_Epmd_index]
max_Epmd_value = maximum(E_p_min_derivative_values) * 1.05 # Adjusted for better annotation placement

# Plotting Fig2b
Fig2b = plot(mu_values, VOI_values, label="Value of Information (VOI)", color=:red, linewidth=2, linestyle=:solid, xlabel="", ylabel="", legend=:bottom, size=(800, 600), framestyle=:origin, legendfontsize=14, tickfontsize=14)
plot!(mu_values, E_p_derivative_values, label="dE[p]/dμ", color=:blue, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_min_derivative_values, label="dE[p_min]/dμ", color=:green, linewidth=2, linestyle=:dashdot)

# Vertical line for max_VOI_mu
max_y_for_line_voi = maximum(VOI_values) # Adjust as necessary
max_y_for_line_Epd = maximum(E_p_derivative_values) # Adjust as necessary
max_y_for_line_Epmd = maximum(E_p_min_derivative_values) # Adjust as necessary
#plot!([max_VOI_mu, max_VOI_mu], [-1.0, max_y_for_line], color=:black, linestyle=:solid, linewidth=2, label="")

# Adjusting annotation to appear on the x-axis
annotate!(Fig2b, max_VOI_mu, max_y_for_line_voi , text(latexstring("$(round(max_VOI_mu, digits=2))"), 14, :center, :bottom))
annotate!(Fig2b, max_Epd_mu, max_y_for_line_Epd , text(latexstring("$(round(max_Epd_mu, digits=2))"), 14, :center, :bottom))
annotate!(Fig2b, max_Epmd_mu, max_y_for_line_Epmd , text(latexstring("$(round(max_Epmd_mu, digits=2))"), 14, :center, :bottom))
annotate!(Fig2b, 0.95, -0.05, text("μ", :right, 14))

# Show plot
display(Fig2b)

# Export to PDF
savefig(Fig2b, join([path_out, "Fig2b.pdf"], "/"))






## Appendix Figure 2, (a), (b), (c), (d) ##

# Parameters for App.Fig2a
N = 2
λ = 0.2
z_max = 0.5
epsilon = 0.03
mu_values = epsilon:0.01:1-epsilon

# Calculating values for Fig1a plots
z_values = [z(mu, λ, z_max) for mu in mu_values]

v = 1.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_1 = E_p_values - E_p_min_values
E_p_derivative_values_1 = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values_1 = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]

v = 2.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_2 = E_p_values - E_p_min_values
E_p_derivative_values_2 = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values_2 = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]

v = 5.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_5 = E_p_values - E_p_min_values
E_p_derivative_values_5 = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values_5 = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]



# Find the μ value that maximizes VOI
max_VOI_index_1 = argmax(VOI_values_1)
max_VOI_mu_1 = mu_values[max_VOI_index_1]
max_VOI_value_1 = maximum(VOI_values_1) * 1.05 # Adjusted for better annotation placement
max_VOI_index_2 = argmax(VOI_values_2)
max_VOI_mu_2 = mu_values[max_VOI_index_2]
max_VOI_value_2 = maximum(VOI_values_2) * 1.05 # Adjusted for better annotation placement
max_VOI_index_5 = argmax(VOI_values_5)
max_VOI_mu_5 = mu_values[max_VOI_index_5]
max_VOI_value_5 = maximum(VOI_values_5) * 1.05 # Adjusted for better annotation placement


# Find the μ value that maximizes E_p_derivative
max_Epd_index_1 = argmax(E_p_derivative_values_1)
max_Epd_mu_1 = mu_values[max_Epd_index_1]
max_Epd_value_1 = maximum(E_p_derivative_values_1) * 1.05 # Adjusted for better annotation placement
max_Epd_index_2 = argmax(E_p_derivative_values_2)
max_Epd_mu_2 = mu_values[max_Epd_index_2]
max_Epd_value_2 = maximum(E_p_derivative_values_2) * 1.05 # Adjusted for better annotation placement
max_Epd_index_5 = argmax(E_p_derivative_values_5)
max_Epd_mu_5 = mu_values[max_Epd_index_5]
max_Epd_value_5 = maximum(E_p_derivative_values_5) * 1.05 # Adjusted for better annotation placement

# Find the μ value that maximizes E_p_min_derivative
max_Epmd_index_1 = argmax(E_p_min_derivative_values_1)
max_Epmd_mu_1 = mu_values[max_Epmd_index_1]
max_Epmd_value_1 = maximum(E_p_min_derivative_values_1) * 1.05 # Adjusted for better annotation placement
max_Epmd_index_2 = argmax(E_p_min_derivative_values_2)
max_Epmd_mu_2 = mu_values[max_Epmd_index_2]
max_Epmd_value_2 = maximum(E_p_min_derivative_values_2) * 1.05 # Adjusted for better annotation placement
max_Epmd_index_5 = argmax(E_p_min_derivative_values_5)
max_Epmd_mu_5 = mu_values[max_Epmd_index_5]
max_Epmd_value_5 = maximum(E_p_min_derivative_values_5) * 1.05 # Adjusted for better annotation placement


# Plotting Appendix Fig2a
AppFig2a = plot(mu_values, VOI_values_1, label="VOI v = 1", color=:red, linewidth=2, linestyle=:solid, xlabel="", ylabel="", legend=:bottom, size=(800, 600), framestyle=:origin, legendfontsize=14, tickfontsize=14)
plot!(mu_values, VOI_values_2, label="VOI v = 2", color=:green, linewidth=2, linestyle=:solid)
plot!(mu_values, VOI_values_5, label="VOI v = 5", color=:purple, linewidth=2, linestyle=:solid)
plot!(mu_values, E_p_derivative_values_1, label="dE[p]/dμ v = 1", color=:red, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_derivative_values_2, label="dE[p]/dμ v = 2", color=:green, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_derivative_values_5, label="dE[p]/dμ v = 5", color=:purple, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_min_derivative_values_1, label="dE[p_min]/dμ v = 1", color=:red, linewidth=2, linestyle=:dot)
plot!(mu_values, E_p_min_derivative_values_2, label="dE[p_min]/dμ v = 2", color=:green, linewidth=2, linestyle=:dot)
plot!(mu_values, E_p_min_derivative_values_5, label="dE[p_min]/dμ v = 5", color=:purple, linewidth=2, linestyle=:dot)

# Vertical line for max_VOI_mu
max_y_for_line_voi_1 = maximum(VOI_values_1) # Adjust as necessary
max_y_for_line_voi_2 = maximum(VOI_values_2) # Adjust as necessary
max_y_for_line_voi_5 = maximum(VOI_values_5) # Adjust as necessary
max_y_for_line_Epd_1 = maximum(E_p_derivative_values_1) # Adjust as necessary
max_y_for_line_Epd_2 = maximum(E_p_derivative_values_2) # Adjust as necessary
max_y_for_line_Epd_5 = maximum(E_p_derivative_values_5) # Adjust as necessary
max_y_for_line_Epmd_1 = maximum(E_p_min_derivative_values_1) # Adjust as necessary
max_y_for_line_Epmd_2 = maximum(E_p_min_derivative_values_2) # Adjust as necessary
max_y_for_line_Epmd_5 = maximum(E_p_min_derivative_values_5) # Adjust as necessary
plot!([max_VOI_mu_1, max_VOI_mu_1], [-4.5, max_y_for_line_voi_5], color=:black, linestyle=:solid, linewidth=2, label="")
plot!([max_Epd_mu_1, max_Epd_mu_1], [-4.5, max_y_for_line_voi_5], color=:black, linestyle=:solid, linewidth=2, label="")
plot!([max_Epmd_mu_1, max_Epmd_mu_1], [-4.5, max_y_for_line_voi_5], color=:black, linestyle=:solid, linewidth=2, label="")

# Adjusting annotation to appear on the x-axis
annotate!(AppFig2a, max_VOI_mu_5, max_y_for_line_voi_5 , text(latexstring("$(round(max_VOI_mu_5, digits=2))"), 14, :center, :bottom))
annotate!(AppFig2a, max_Epd_mu_1, max_y_for_line_voi_5 , text(latexstring("$(round(max_Epd_mu_5, digits=2))"), 14, :center, :bottom))
annotate!(AppFig2a, max_Epmd_mu_1, max_y_for_line_voi_5 , text(latexstring("$(round(max_Epmd_mu_5, digits=2))"), 14, :center, :bottom))
annotate!(AppFig2a, 0.95, -0.28, text("μ", :right, 14))

# Show plot
display(AppFig2a)

# Export to PDF
savefig(AppFig2a, join([path_out, "AppFig2a.pdf"], "/"))




# Parameters for App.Fig2b
N = 3
λ = 0.2
z_max = 0.5
epsilon = 0.04
mu_values = epsilon:0.01:1-epsilon

# Calculating values for Fig1a plots
z_values = [z(mu, λ, z_max) for mu in mu_values]

v = 1.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_1 = E_p_values - E_p_min_values
E_p_derivative_values_1 = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values_1 = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]

v = 2.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_2 = E_p_values - E_p_min_values
E_p_derivative_values_2 = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values_2 = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]

v = 5.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_5 = E_p_values - E_p_min_values
E_p_derivative_values_5 = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values_5 = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]



# Find the μ value that maximizes VOI
max_VOI_index_1 = argmax(VOI_values_1)
max_VOI_mu_1 = mu_values[max_VOI_index_1]
max_VOI_value_1 = maximum(VOI_values_1) * 1.05 # Adjusted for better annotation placement
max_VOI_index_2 = argmax(VOI_values_2)
max_VOI_mu_2 = mu_values[max_VOI_index_2]
max_VOI_value_2 = maximum(VOI_values_2) * 1.05 # Adjusted for better annotation placement
max_VOI_index_5 = argmax(VOI_values_5)
max_VOI_mu_5 = mu_values[max_VOI_index_5]
max_VOI_value_5 = maximum(VOI_values_5) * 1.05 # Adjusted for better annotation placement


# Find the μ value that maximizes E_p_derivative
max_Epd_index_1 = argmax(E_p_derivative_values_1)
max_Epd_mu_1 = mu_values[max_Epd_index_1]
max_Epd_value_1 = maximum(E_p_derivative_values_1) * 1.05 # Adjusted for better annotation placement
max_Epd_index_2 = argmax(E_p_derivative_values_2)
max_Epd_mu_2 = mu_values[max_Epd_index_2]
max_Epd_value_2 = maximum(E_p_derivative_values_2) * 1.05 # Adjusted for better annotation placement
max_Epd_index_5 = argmax(E_p_derivative_values_5)
max_Epd_mu_5 = mu_values[max_Epd_index_5]
max_Epd_value_5 = maximum(E_p_derivative_values_5) * 1.05 # Adjusted for better annotation placement

# Find the μ value that maximizes E_p_min_derivative
max_Epmd_index_1 = argmax(E_p_min_derivative_values_1)
max_Epmd_mu_1 = mu_values[max_Epmd_index_1]
max_Epmd_value_1 = maximum(E_p_min_derivative_values_1) * 1.05 # Adjusted for better annotation placement
max_Epmd_index_2 = argmax(E_p_min_derivative_values_2)
max_Epmd_mu_2 = mu_values[max_Epmd_index_2]
max_Epmd_value_2 = maximum(E_p_min_derivative_values_2) * 1.05 # Adjusted for better annotation placement
max_Epmd_index_5 = argmax(E_p_min_derivative_values_5)
max_Epmd_mu_5 = mu_values[max_Epmd_index_5]
max_Epmd_value_5 = maximum(E_p_min_derivative_values_5) * 1.05 # Adjusted for better annotation placement


# Plotting Appendix Fig2b
AppFig2b = plot(mu_values, VOI_values_1, label="VOI v = 1", color=:red, linewidth=2, linestyle=:solid, xlabel="", ylabel="", legend=:bottom, size=(800, 600), framestyle=:origin, legendfontsize=14, tickfontsize=14)
plot!(mu_values, VOI_values_2, label="VOI v = 2", color=:green, linewidth=2, linestyle=:solid)
plot!(mu_values, VOI_values_5, label="VOI v = 5", color=:purple, linewidth=2, linestyle=:solid)
plot!(mu_values, E_p_derivative_values_1, label="dE[p]/dμ v = 1", color=:red, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_derivative_values_2, label="dE[p]/dμ v = 2", color=:green, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_derivative_values_5, label="dE[p]/dμ v = 5", color=:purple, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_min_derivative_values_1, label="dE[p_min]/dμ v = 1", color=:red, linewidth=2, linestyle=:dot)
plot!(mu_values, E_p_min_derivative_values_2, label="dE[p_min]/dμ v = 2", color=:green, linewidth=2, linestyle=:dot)
plot!(mu_values, E_p_min_derivative_values_5, label="dE[p_min]/dμ v = 5", color=:purple, linewidth=2, linestyle=:dot)

# Vertical line for max_VOI_mu
max_y_for_line_voi_1 = maximum(VOI_values_1) # Adjust as necessary
max_y_for_line_voi_2 = maximum(VOI_values_2) # Adjust as necessary
max_y_for_line_voi_5 = maximum(VOI_values_5) # Adjust as necessary
max_y_for_line_Epd_1 = maximum(E_p_derivative_values_1) # Adjust as necessary
max_y_for_line_Epd_2 = maximum(E_p_derivative_values_2) # Adjust as necessary
max_y_for_line_Epd_5 = maximum(E_p_derivative_values_5) # Adjust as necessary
max_y_for_line_Epmd_1 = maximum(E_p_min_derivative_values_1) # Adjust as necessary
max_y_for_line_Epmd_2 = maximum(E_p_min_derivative_values_2) # Adjust as necessary
max_y_for_line_Epmd_5 = maximum(E_p_min_derivative_values_5) # Adjust as necessary
plot!([max_VOI_mu_1, max_VOI_mu_1], [-5.2, max_y_for_line_voi_5], color=:black, linestyle=:solid, linewidth=2, label="")
plot!([max_Epd_mu_1, max_Epd_mu_1], [-5.2, max_y_for_line_voi_5], color=:black, linestyle=:solid, linewidth=2, label="")
plot!([max_Epmd_mu_1, max_Epmd_mu_1], [-5.2, max_y_for_line_voi_5], color=:black, linestyle=:solid, linewidth=2, label="")

# Adjusting annotation to appear on the x-axis
annotate!(AppFig2b, max_VOI_mu_5, max_y_for_line_voi_5 , text(latexstring("$(round(max_VOI_mu_5, digits=2))"), 14, :center, :bottom))
annotate!(AppFig2b, max_Epd_mu_1, max_y_for_line_voi_5 , text(latexstring("$(round(max_Epd_mu_5, digits=2))"), 14, :center, :bottom))
annotate!(AppFig2b, max_Epmd_mu_1, max_y_for_line_voi_5 , text(latexstring("$(round(max_Epmd_mu_5, digits=2))"), 14, :center, :bottom))
annotate!(AppFig2b, 0.95, -0.3, text("μ", :right, 14))

# Show plot
display(AppFig2b)

# Export to PDF
savefig(AppFig2b, join([path_out, "AppFig2b.pdf"], "/"))



# Parameters for App.Fig2c
N = 5
λ = 0.2
z_max = 0.5
epsilon = 0.05
mu_values = epsilon:0.01:1-epsilon

# Calculating values for Fig2c plots
z_values = [z(mu, λ, z_max) for mu in mu_values]

v = 1.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_1 = E_p_values - E_p_min_values
E_p_derivative_values_1 = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values_1 = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]

v = 2.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_2 = E_p_values - E_p_min_values
E_p_derivative_values_2 = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values_2 = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]

v = 5.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_5 = E_p_values - E_p_min_values
E_p_derivative_values_5 = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values_5 = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]



# Find the μ value that maximizes VOI
max_VOI_index_1 = argmax(VOI_values_1)
max_VOI_mu_1 = mu_values[max_VOI_index_1]
max_VOI_value_1 = maximum(VOI_values_1) * 1.05 # Adjusted for better annotation placement
max_VOI_index_2 = argmax(VOI_values_2)
max_VOI_mu_2 = mu_values[max_VOI_index_2]
max_VOI_value_2 = maximum(VOI_values_2) * 1.05 # Adjusted for better annotation placement
max_VOI_index_5 = argmax(VOI_values_5)
max_VOI_mu_5 = mu_values[max_VOI_index_5]
max_VOI_value_5 = maximum(VOI_values_5) * 1.05 # Adjusted for better annotation placement


# Find the μ value that maximizes E_p_derivative
max_Epd_index_1 = argmax(E_p_derivative_values_1)
max_Epd_mu_1 = mu_values[max_Epd_index_1]
max_Epd_value_1 = maximum(E_p_derivative_values_1) * 1.05 # Adjusted for better annotation placement
max_Epd_index_2 = argmax(E_p_derivative_values_2)
max_Epd_mu_2 = mu_values[max_Epd_index_2]
max_Epd_value_2 = maximum(E_p_derivative_values_2) * 1.05 # Adjusted for better annotation placement
max_Epd_index_5 = argmax(E_p_derivative_values_5)
max_Epd_mu_5 = mu_values[max_Epd_index_5]
max_Epd_value_5 = maximum(E_p_derivative_values_5) * 1.05 # Adjusted for better annotation placement

# Find the μ value that maximizes E_p_min_derivative
max_Epmd_index_1 = argmax(E_p_min_derivative_values_1)
max_Epmd_mu_1 = mu_values[max_Epmd_index_1]
max_Epmd_value_1 = maximum(E_p_min_derivative_values_1) * 1.05 # Adjusted for better annotation placement
max_Epmd_index_2 = argmax(E_p_min_derivative_values_2)
max_Epmd_mu_2 = mu_values[max_Epmd_index_2]
max_Epmd_value_2 = maximum(E_p_min_derivative_values_2) * 1.05 # Adjusted for better annotation placement
max_Epmd_index_5 = argmax(E_p_min_derivative_values_5)
max_Epmd_mu_5 = mu_values[max_Epmd_index_5]
max_Epmd_value_5 = maximum(E_p_min_derivative_values_5) * 1.05 # Adjusted for better annotation placement


# Plotting Appendix Fig2c
AppFig2c = plot(mu_values, VOI_values_1, label="VOI v = 1", color=:red, linewidth=2, linestyle=:solid, xlabel="", ylabel="", legend=:bottom, size=(800, 600), framestyle=:origin, legendfontsize=14, tickfontsize=14)
plot!(mu_values, VOI_values_2, label="VOI v = 2", color=:green, linewidth=2, linestyle=:solid)
plot!(mu_values, VOI_values_5, label="VOI v = 5", color=:purple, linewidth=2, linestyle=:solid)
plot!(mu_values, E_p_derivative_values_1, label="dE[p]/dμ v = 1", color=:red, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_derivative_values_2, label="dE[p]/dμ v = 2", color=:green, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_derivative_values_5, label="dE[p]/dμ v = 5", color=:purple, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_min_derivative_values_1, label="dE[p_min]/dμ v = 1", color=:red, linewidth=2, linestyle=:dot)
plot!(mu_values, E_p_min_derivative_values_2, label="dE[p_min]/dμ v = 2", color=:green, linewidth=2, linestyle=:dot)
plot!(mu_values, E_p_min_derivative_values_5, label="dE[p_min]/dμ v = 5", color=:purple, linewidth=2, linestyle=:dot)

# Vertical line for max_VOI_mu
max_y_for_line_voi_1 = maximum(VOI_values_1) # Adjust as necessary
max_y_for_line_voi_2 = maximum(VOI_values_2) # Adjust as necessary
max_y_for_line_voi_5 = maximum(VOI_values_5) # Adjust as necessary
max_y_for_line_Epd_1 = maximum(E_p_derivative_values_1) # Adjust as necessary
max_y_for_line_Epd_2 = maximum(E_p_derivative_values_2) # Adjust as necessary
max_y_for_line_Epd_5 = maximum(E_p_derivative_values_5) # Adjust as necessary
max_y_for_line_Epmd_1 = maximum(E_p_min_derivative_values_1) # Adjust as necessary
max_y_for_line_Epmd_2 = maximum(E_p_min_derivative_values_2) # Adjust as necessary
max_y_for_line_Epmd_5 = maximum(E_p_min_derivative_values_5) # Adjust as necessary
plot!([max_VOI_mu_1, max_VOI_mu_1], [-5.25, max_y_for_line_voi_5], color=:black, linestyle=:solid, linewidth=2, label="")
plot!([max_Epd_mu_1, max_Epd_mu_1], [-5.25, max_y_for_line_voi_5], color=:black, linestyle=:solid, linewidth=2, label="")
plot!([max_Epmd_mu_1, max_Epmd_mu_1], [-5.25, max_y_for_line_voi_5], color=:black, linestyle=:solid, linewidth=2, label="")

# Adjusting annotation to appear on the x-axis
annotate!(AppFig2c, max_VOI_mu_5, max_y_for_line_voi_5 , text(latexstring("$(round(max_VOI_mu_5, digits=2))"), 14, :center, :bottom))
annotate!(AppFig2c, max_Epd_mu_1, max_y_for_line_voi_5 , text(latexstring("$(round(max_Epd_mu_5, digits=2))"), 14, :center, :bottom))
annotate!(AppFig2c, max_Epmd_mu_1, max_y_for_line_voi_5 , text(latexstring("$(round(max_Epmd_mu_5, digits=2))"), 14, :center, :bottom))
annotate!(AppFig2c, 0.95, -0.4, text("μ", :right, 14))

# Show plot
display(AppFig2c)

# Export to PDF
savefig(AppFig2c, join([path_out, "AppFig2c.pdf"], "/"))




# Parameters for App.Fig2d
N = 10
λ = 0.2
z_max = 0.5
epsilon = 0.05
mu_values = epsilon:0.01:1-epsilon

# Calculating values for Fig2d plots
z_values = [z(mu, λ, z_max) for mu in mu_values]

v = 1.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_1 = E_p_values - E_p_min_values
E_p_derivative_values_1 = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values_1 = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]

v = 2.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_2 = E_p_values - E_p_min_values
E_p_derivative_values_2 = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values_2 = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]

v = 5.0
E_p_values = [E_p(mu, N, v) for mu in mu_values]
E_p_min_values = [E_p_min(mu, E_p_values[i]) for (i, mu) in enumerate(mu_values)]
VOI_values_5 = E_p_values - E_p_min_values
E_p_derivative_values_5 = [numerical_derivative(mu -> E_p(mu, N, v), mu) for mu in mu_values]
E_p_min_derivative_values_5 = [numerical_derivative(mu -> E_p_min(mu, E_p(mu, N, v)), mu) for mu in mu_values]



# Find the μ value that maximizes VOI
max_VOI_index_1 = argmax(VOI_values_1)
max_VOI_mu_1 = mu_values[max_VOI_index_1]
max_VOI_value_1 = maximum(VOI_values_1) * 1.05 # Adjusted for better annotation placement
max_VOI_index_2 = argmax(VOI_values_2)
max_VOI_mu_2 = mu_values[max_VOI_index_2]
max_VOI_value_2 = maximum(VOI_values_2) * 1.05 # Adjusted for better annotation placement
max_VOI_index_5 = argmax(VOI_values_5)
max_VOI_mu_5 = mu_values[max_VOI_index_5]
max_VOI_value_5 = maximum(VOI_values_5) * 1.05 # Adjusted for better annotation placement


# Find the μ value that maximizes E_p_derivative
max_Epd_index_1 = argmax(E_p_derivative_values_1)
max_Epd_mu_1 = mu_values[max_Epd_index_1]
max_Epd_value_1 = maximum(E_p_derivative_values_1) * 1.05 # Adjusted for better annotation placement
max_Epd_index_2 = argmax(E_p_derivative_values_2)
max_Epd_mu_2 = mu_values[max_Epd_index_2]
max_Epd_value_2 = maximum(E_p_derivative_values_2) * 1.05 # Adjusted for better annotation placement
max_Epd_index_5 = argmax(E_p_derivative_values_5)
max_Epd_mu_5 = mu_values[max_Epd_index_5]
max_Epd_value_5 = maximum(E_p_derivative_values_5) * 1.05 # Adjusted for better annotation placement

# Find the μ value that maximizes E_p_min_derivative
max_Epmd_index_1 = argmax(E_p_min_derivative_values_1)
max_Epmd_mu_1 = mu_values[max_Epmd_index_1]
max_Epmd_value_1 = maximum(E_p_min_derivative_values_1) * 1.05 # Adjusted for better annotation placement
max_Epmd_index_2 = argmax(E_p_min_derivative_values_2)
max_Epmd_mu_2 = mu_values[max_Epmd_index_2]
max_Epmd_value_2 = maximum(E_p_min_derivative_values_2) * 1.05 # Adjusted for better annotation placement
max_Epmd_index_5 = argmax(E_p_min_derivative_values_5)
max_Epmd_mu_5 = mu_values[max_Epmd_index_5]
max_Epmd_value_5 = maximum(E_p_min_derivative_values_5) * 1.05 # Adjusted for better annotation placement


# Plotting Appendix Fig2d
AppFig2d = plot(mu_values, VOI_values_1, label="VOI v = 1", color=:red, linewidth=2, linestyle=:solid, xlabel="", ylabel="", legend=:bottom, size=(800, 600), framestyle=:origin, legendfontsize=14, tickfontsize=14)
plot!(mu_values, VOI_values_2, label="VOI v = 2", color=:green, linewidth=2, linestyle=:solid)
plot!(mu_values, VOI_values_5, label="VOI v = 5", color=:purple, linewidth=2, linestyle=:solid)
plot!(mu_values, E_p_derivative_values_1, label="dE[p]/dμ v = 1", color=:red, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_derivative_values_2, label="dE[p]/dμ v = 2", color=:green, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_derivative_values_5, label="dE[p]/dμ v = 5", color=:purple, linewidth=2, linestyle=:dash)
plot!(mu_values, E_p_min_derivative_values_1, label="dE[p_min]/dμ v = 1", color=:red, linewidth=2, linestyle=:dot)
plot!(mu_values, E_p_min_derivative_values_2, label="dE[p_min]/dμ v = 2", color=:green, linewidth=2, linestyle=:dot)
plot!(mu_values, E_p_min_derivative_values_5, label="dE[p_min]/dμ v = 5", color=:purple, linewidth=2, linestyle=:dot)

# Vertical line for max_VOI_mu
max_y_for_line_voi_1 = maximum(VOI_values_1) # Adjust as necessary
max_y_for_line_voi_2 = maximum(VOI_values_2) # Adjust as necessary
max_y_for_line_voi_5 = maximum(VOI_values_5) # Adjust as necessary
max_y_for_line_Epd_1 = maximum(E_p_derivative_values_1) # Adjust as necessary
max_y_for_line_Epd_2 = maximum(E_p_derivative_values_2) # Adjust as necessary
max_y_for_line_Epd_5 = maximum(E_p_derivative_values_5) # Adjust as necessary
max_y_for_line_Epmd_1 = maximum(E_p_min_derivative_values_1) # Adjust as necessary
max_y_for_line_Epmd_2 = maximum(E_p_min_derivative_values_2) # Adjust as necessary
max_y_for_line_Epmd_5 = maximum(E_p_min_derivative_values_5) # Adjust as necessary
plot!([max_VOI_mu_1, max_VOI_mu_1], [-7.7, max_y_for_line_voi_5], color=:black, linestyle=:solid, linewidth=2, label="")
plot!([max_Epd_mu_1, max_Epd_mu_1], [-7.7, max_y_for_line_voi_5], color=:black, linestyle=:solid, linewidth=2, label="")
plot!([max_Epmd_mu_1, max_Epmd_mu_1], [-7.7, max_y_for_line_voi_5], color=:black, linestyle=:solid, linewidth=2, label="")

# Adjusting annotation to appear on the x-axis
annotate!(AppFig2d, max_VOI_mu_5, max_y_for_line_voi_5 , text(latexstring("$(round(max_VOI_mu_5, digits=2))"), 14, :center, :bottom))
annotate!(AppFig2d, max_Epd_mu_1, max_y_for_line_voi_5 , text(latexstring("$(round(max_Epd_mu_5, digits=2))"), 14, :center, :bottom))
annotate!(AppFig2d, max_Epmd_mu_1, max_y_for_line_voi_5 , text(latexstring("$(round(max_Epmd_mu_5, digits=2))"), 14, :center, :bottom))
annotate!(AppFig2d, 0.95, -0.65, text("μ", :right, 14))

# Show plot
display(AppFig2d)

# Export to PDF
savefig(AppFig2d, join([path_out, "AppFig2d.pdf"], "/"))